# Author: Xinwu
# Describe: 
# Completion Time: 
# Email: lexinwu@outlook.com

library(Seurat)
library(monocle)
source("function.R")
######################################################### 标准分析 ########################################################
# get data from project.rdata
result <- prep_data(project = project)

# create monocle object
monocle_obj <- newCellDataSet(cellData = methods::as(as.matrix(result$count), "sparseMatrix"),
                              phenoData = new("AnnotatedDataFrame", data = result$pheno),
                              featureData = new("AnnotatedDataFrame", data = result$anno),
                              lowerDetectionLimit = 0.1,
                              expressionFamily = negbinomial.size())

# 计算(评估) SizeFactor 和 Dispersion(离散度)
monocle_obj <- estimateSizeFactors(monocle_obj)
monocle_obj <- estimateDispersions(monocle_obj)

# 过滤细胞(可选, 如果是已经跑过了 Seurat 标准流程, 就不需要)

# 选择特征基因, 轨迹推断的第一步就是选择 monocle 将用作机器学习的基因, 也叫特征选择, 它对轨迹的形状有很大影响
# 1. 自定义发育 marker 基因, 从文件读入即可。
# 2. 选择差异基因(cluster、celltype)

# 3. 选择高变基因(Seurat、monocle)
# Seurat 计算高变基因
project <- FindVariableFeatures(project, selection.method = 'vst', nfeatures = 2000)
ordergene <- VariableFeatures(project)
# monocle 计算高变基因
monocle_obj <- estimateDispersions(monocle_obj)
disp_table <- dispersionTable(monocle_obj)
ordergene <- subset(disp_table,mean_expression >= 0.3 &  dispersion_empirical >= 1 * dispersion_fit)$gene_id
monocle_obj <- setOrderingFilter(monocle_obj, ordergene)

# 降维
monocle_obj <- reduceDimension(monocle_obj, max_components = 2, method = 'DDRTree', norm_method = c("log"))

# 轨迹构建 和 重设 root state
monocle_obj <- orderCells(monocle_obj)
monocle_obj <- orderCells(monocle_obj, root_stat = 2) # 自定义的 root state

# 轨迹可视化, monocle 使用的是 ggplot2 来进行绘图的, 所以可以使用 ggplot2 的方式进行图片的修改(分面、颜色、主题 等都可以自定义修改)
# 可以使用 colnames(monocle_obj@phenoData@data) 里的值对轨迹进行着色
# fData(monocle_obj) 获取 monocle_obj 对象中的 基因描述信息, 同 monocle_obj@featureData@data
# pData(monocle_obj) 获取 monocle_obj 对象中的 细胞描述信息, 同 monocle_obj@phenoData@data
plot_cell_trajectory(monocle_obj, color_by = , )

# 基因在伪时间上的表达变化情况
gene_list <- c()
plot_genes_in_pseudotime(monocl_obj[gene_list,], color_by = ,)

# 在轨迹图上展示单个基因的表达量 如 CCl5
pData(monocle_obj)$CCl5 <- log2(exprs(monocl_obj)['CCl5',] + 1) # 把 CCl5 的表达量信息, 附加在 monocle_obj@phenoData@data 中
plot_cell_trajectory(monocle_obj, color_by = "CCl5")

# 差异基因 sm.ns 函数指出 monocle 应该通过表达值拟合自然样条曲线, 以帮助它将表达式的变化描述为进程的函数, 好像只要在接 Pseudotime 的时候会使用
diff <- differentialGeneTest(cds = monocl_obj, fullModelFormulaStr = "~sm.ns(Pseudotime)", cores = 4) # 拟时差异基因
diff <- differentialGeneTest(cds = monocl_obj, fullModelFormulaStr = "~State", cores = 4) # State差异基因

# 挑选部分差异基因进行热图展示
sig_gene_names <- rownames(diff)[1:100]
plot_pseudotime_heatmap(monocl_obj[sig_gene_names,], num_clusters = 3, cores = 1, show_rownames = T,return_heatmap = T)

# 分支分析(BEAM)
BEAM_res <- BEAM(monocl_obj, branch_point = 1, cores = 1) # branch_point 就是轨迹图上 黑色圆圈上的数字
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
plot_genes_branched_heatmap(monocl_obj[row.names(subset(BEAM_res,
                                          qval < 1e-4)),],
                                          branch_point = 1,
                                          num_clusters = 4,
                                          cores = 1,
                                          use_gene_short_name = T,
                                          show_rownames = T)

# 选择分支分析中感兴趣的基因进行可视化, 注意进行单个基因展示时会报错
genes <- row.names(subset(fData(monocl_obj),
          gene_short_name %in% c("Ccnd2", "Sftpb", "Pdpn")))
plot_genes_branched_pseudotime(monocl_obj[genes,],
                       branch_point = 1,
                       color_by = "State",
                       ncol = 1)

# 个性化可视化
# 1. 轨迹图上加饼图
# p1 为以 Cluster 进行着色的轨迹图, p2~p4 为 state1~3 细胞占比的饼图, 然后把这4张图进行 merge 即可
p1 <- plot_cell_trajectory(monocl_data,color_by="Cluster") + 
  scale_color_manual(breaks = c(2,17), 
                     values=c("#f4c768","#79bfdb"),
                     guide = ggplot2::guide_legend(override.aes = list(size = 6))) + 
  blank_theme +
  theme(legend.position = "top") 
p2 <- plot_pie(monocl_obj = monocl_data,filter_by = "State",
               retain_by = "1",group_by = "Cluster",
               col = c("#f4c768","#79bfdb"))+
  theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"))
p3 <- plot_pie(monocl_obj = monocl_data,filter_by = "State",
               retain_by = "2",group_by = "Cluster",
               col = c("#f4c768","#79bfdb"))+
  theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"))
p4 <- plot_pie(monocl_obj = monocl_data,filter_by = "State",
               retain_by = "3",group_by = "Cluster",
               col = c("#f4c768","#79bfdb"))+
  theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"))

p <- merge_fig(main_fig = p1, anno_fig = p2, x = -2, y = -1, d = 2)
p <- merge_fig(main_fig = p, anno_fig = p3, x = 3, y = 1.5, d = 2)
p <- merge_fig(main_fig = p, anno_fig = p4, x = 2.2, y = -1.5, d = 2)